Certificado de Aceptacion:
Study to Determine Levels of Cadmiumin Cocoa Crops Applied to Inland Areas of Peru
Paper:
Acreditation
El link es el siguiente: Agronomy MDPI.
El link es el siguiente: ECITEC UNI 2020
DEM <- raster("Data_Suelos/raster/DEM_DRONE/Zona_Oeste.tif")
# str(DEM)
# DEM
# slope <- terrain(DEM, opt ="slope")
# aspect <- terrain(DEM, opt="aspect")
# DEM.hill <- hillShade(slope, aspect, angle = 40, direction = 270)
DEM <- crop(DEM, extent(518600, 519200, 9031800, 9032400))
#plot(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
# zlim=c(0, 250))
image(DEM, main = "Drone - Digital Elevation Model (DEM)", col = topo.colors(20),
zlim=c(0, 250))## Reading layer `Age_Planta' from data source
## `D:\CursoDAG\DataAnalisisGeociencias\Modulos\ModuloII\Data_Suelos\shp\Cacao_Factors\Age_Planta.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 6 features and 5 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 518541 ymin: 9031841 xmax: 519584 ymax: 9032609
## Projected CRS: WGS 84 / UTM zone 18S
mapview(DEM, legend =TRUE)+
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10,
lwd = 2, color = "blue")## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ;
## the supplied Raster* has 573840800
## ... decreasing Raster* resolution to 5e+05 pixels
## to view full resolution set 'maxpixels = 573840800 '
Factores_cacao2 <- as.data.frame(Factores_cacao)
#Factores_cacao2$Uso_Tierra <- edit(Factores_cacao2$Uso_Tierra)
Factores_cacao2$Var_cocoa <- c("Creole-Aromatic","CCN51","","CCN51","Creole-Aromatic",
"CCN51")
Factores_cacao3 <- st_as_sf(Factores_cacao2)
mapview(Factores_cacao3, zcol="Var_cocoa")Honoria_elementos <- readxl::read_xlsx(path = "Data_Suelos/csv_xlsx_txt/BD_Cd.xlsx")
str(Honoria_elementos)## tibble [24 × 7] (S3: tbl_df/tbl/data.frame)
## $ Xeast : num [1:24] 518642 518644 518690 518701 518706 ...
## $ Ynorth : num [1:24] 9032291 9032343 9032334 9032232 9032284 ...
## $ Cdsoil_mgkg : num [1:24] 1.75 1.6 1.9 1.7 1.25 1.75 1.5 1.4 1.4 1.35 ...
## $ Pbsoil_mgkg : num [1:24] 5.41 5.03 7.24 5.81 5.91 5.62 6.57 5.26 5.01 7.02 ...
## $ Znsoil_mgkg : num [1:24] 70.5 99.4 77.1 97.1 74.7 ...
## $ pHapprox : num [1:24] 7.06 6.55 6.28 8 6.7 6.4 6.53 5.79 5.41 6.39 ...
## $ Ceapprox_uS_cm: num [1:24] 39.9 10.3 5 126.7 8 ...
Honoria_elementos <- Honoria_elementos %>% rename(Este= Xeast, Norte=Ynorth,
Cd = Cdsoil_mgkg,
Pb = Pbsoil_mgkg,
Zn = Znsoil_mgkg,
pH = pHapprox,
Ce = Ceapprox_uS_cm)
Honoria_elementos2 <-Honoria_elementos[ ,c("Norte","Este")]
Honoria_elementos2 <- Honoria_elementos2[ ,order(c(names(Honoria_elementos2)))]
sputm <- SpatialPoints(Honoria_elementos2, proj4string=CRS("+proj=utm +zone=18 +south +datum=WGS84"))
spgeo <- spTransform(sputm, CRS("+proj=longlat +datum=WGS84"))
spgeo <- as.data.frame(spgeo)
colnames(spgeo) <- c("lng","lat")
Honoria_elementos2 <- cbind(Honoria_elementos,spgeo)
Honoria_elementos3 <- st_as_sf(Honoria_elementos2, coords = c("lng", "lat"),
remove = FALSE, crs = 4326, agr = "constant")
mapview(Honoria_elementos3, zcol = "Cd") +mapview(Factores_cacao3, zcol="Var_cocoa")+
mapview(Factores_cacao3, zcol = "Uso_Tierra")+
mapview(Factores_cacao3, zcol = "Age")Factores_cacao4 <- Factores_cacao3 %>%
st_transform(crs = 4326 )
Honoria_elementos3 %>%
st_intersection(Factores_cacao4) %>%
mapview(zcol = "Uso_Tierra", burst = TRUE, cex = "Cd")+
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10,
lwd = 2, color = "blue")Honoria_elementos3 %>%
st_intersection(Factores_cacao4) %>%
mapview(zcol = "Age", burst = TRUE, cex = "Cd" )+
mapview(Factores_cacao, alpha.regions = 0.01 , legend.opacity = 0.10,
lwd = 2, color = "blue")Factores_cacao_n <- Factores_cacao
Honoria_elementos_n <- Honoria_elementos3 %>% st_transform(crs = st_crs(Factores_cacao))
# sf_use_s2(TRUE) #cambiar a spherical geometry
intersection <- st_intersection(Honoria_elementos_n, Factores_cacao_n)
intersectionmin_max_mean_sd <- list(
min = ~min(.x, na.rm = TRUE),
max = ~max(.x, na.rm = TRUE),
mean = ~mean(.x, na.rm = TRUE),
sd = ~sd(.x, na.rm = TRUE)
)Ver vignette("colwise") para más detalles.
correlacion <- cor(Honoria_statical[ ,c("Cd", "Pb", "Zn", "pH", "Ce")]
, method = "pearson")
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria",
hclust.method = "median", addrect = 2)
corrplot(correlacion, method ="circle", diag = TRUE, title= "Correlacion de Pearson Honoria")correlacion2 <- cor(Honoria_statical[ Honoria_statical$Uso_Tierra=="Maiz 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
# correlacion2 <- cor(Honoria_statical %>% filter(Uso_Tierra=="Maiz 10") %>%
# select(Cd, Pb, Zn, pH, Ce), method = "pearson")
corrplot(correlacion2, method ="number")correlacion3 <- cor(Honoria_statical[Honoria_statical$Uso_Tierra=="Pasto 10", c("Cd", "Pb", "Zn", "pH", "Ce")], method = "pearson")
corrplot(correlacion3, method ="circle")Para ver toda la funcionalidad revisar : An Introduction to corrplot Package.
## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa
prcomp()res.pca <- prcomp(Honoria_pca, scale = TRUE)
res.pca <- PCA(Honoria_pca, graph=FALSE)
print(res.pca)## **Results for the Principal Component Analysis (PCA)**
## The analysis was performed on 22 individuals, described by 5 variables
## *The results are available in the following objects:
##
## name description
## 1 "$eig" "eigenvalues"
## 2 "$var" "results for the variables"
## 3 "$var$coord" "coord. for the variables"
## 4 "$var$cor" "correlations variables - dimensions"
## 5 "$var$cos2" "cos2 for the variables"
## 6 "$var$contrib" "contributions of the variables"
## 7 "$ind" "results for the individuals"
## 8 "$ind$coord" "coord. for the individuals"
## 9 "$ind$cos2" "cos2 for the individuals"
## 10 "$ind$contrib" "contributions of the individuals"
## 11 "$call" "summary statistics"
## 12 "$call$centre" "mean of the variables"
## 13 "$call$ecart.type" "standard error of the variables"
## 14 "$call$row.w" "weights for the individuals"
## 15 "$call$col.w" "weights for the variables"
## [1] 1.664091 5.980000 59.403182 6.530000 27.547727
# La varianza explicada por el primer eigenvalor es 38.60% y el segundo 25.58% los cuales
# seran tomados para el analisis de los componentes principales representando 64.19% de la
# varianza total. Según (Kaiser 1961) un eigenvalue>1 indica un punto de corte
# para contener el deja a voluntad del investigador realizar el corte en un valor adecuado
# para la visualizacion.
eig.val<- get_eigenvalue(res.pca)
eig.val## eigenvalue variance.percent cumulative.variance.percent
## Dim.1 1.9300862 38.601723 38.60172
## Dim.2 1.2794547 25.589095 64.19082
## Dim.3 0.9204006 18.408013 82.59883
## Dim.4 0.6054406 12.108812 94.70764
## Dim.5 0.2646179 5.292357 100.00000
## Principal Component Analysis Results for variables
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the variables"
## 2 "$cor" "Correlations between variables and dimensions"
## 3 "$cos2" "Cos2 for the variables"
## 4 "$contrib" "contributions of the variables"
#Circulo de Correlacion:
#La correlacion entre las variables de los componentes principales son usadas como
# coordenadas de la variable para los PC. Segun esto las observaciones son representadas
# por su proyeccion, pero las variables son representadas por sus correlaciones.
var$coord## Dim.1 Dim.2 Dim.3 Dim.4 Dim.5
## Cd 0.5097511 -0.5694786 0.50836826 -0.35210702 0.18283969
## Pb 0.6597640 -0.4324834 0.04309302 0.61090780 -0.05103086
## Zn -0.2583970 0.5477468 0.75441302 0.24677884 0.05626362
## pH 0.7076697 0.5604771 -0.27468642 0.02387721 0.33022158
## Ce 0.8169370 0.3923593 0.12455405 -0.21629371 -0.34113264
# El plot de correlacion nos indica que las variables pH y Ce estan correlacionadas
# positivamente (agrupadas), asi como el Cd, Pb. Las variables Cd-Zn y Pb-Zn tienen
# una correlación inversa, mientras que entre ph-Zn y Ce-Zn no existe una relacion
# significativa.
# Las variables pH, Ce, Pb y Cd estan bien representadas, adicionalmente
# se debe notar que el Zn disminuye su representatividad.##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## Ce 0.8169370 3.498136e-06
## pH 0.7076697 2.296007e-04
## Pb 0.6597640 8.356745e-04
## Cd 0.5097511 1.537322e-02
##
## Link between the variable and the continuous variables (R-square)
## =================================================================================
## correlation p.value
## pH 0.5604771 0.006664419
## Zn 0.5477468 0.008319053
## Pb -0.4324834 0.044405255
## Cd -0.5694786 0.005667244
## Principal Component Analysis Results for individuals
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the individuals"
## 2 "$cos2" "Cos2 for the individuals"
## 3 "$contrib" "contributions of the individuals"
fviz_pca_ind(res.pca,
geom.ind = "point",
col.ind = Honoria_pca_var$Age,
palette = rainbow(6),
addEllipses = TRUE,
legend.title="Age")## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
fviz_pca_ind(res.pca,
geom.ind = "point",
col.ind = Honoria_pca_var$Uso_Tierra,
palette = rainbow(4),
addEllipses = TRUE,
legend.title="Land Use")## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
q <- fviz_pca_biplot(res.pca,
geom.ind ="point",
fill.ind = Honoria_pca_var$Age, col.ind = "black",
pointshape = 21, pointsize = 2,arrowsize=1,
palette = rainbow(6),
addEllipses = TRUE,
ellipse.level= 0.95,
col.var = "black",
gradient.cols = "RdYlBu",
legend.title = "Time of \n Age",
repel=FALSE
)
q## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse
# Compute PCA with ncp = 3
res.pca1 <- PCA(Honoria_pca, ncp = 6, graph = FALSE)
# Compute hierarchical clustering on principal components
res.hcpc <- HCPC(res.pca1, graph = FALSE)fviz_dend(res.hcpc,
cex = 0.7, # Label size
palette = "jco", # Color palette see ?ggpubr::ggpar
rect = TRUE, rect_fill = TRUE, # Add rectangle around groups
rect_border = "jco", # Rectangle color
labels_track_height = 0.8 # Augment the room for labels
)## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <]8;;https://github.com/kassambara/factoextra/issueshttps://github.com/kassambara/factoextra/issues]8;;>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
fviz_cluster(res.hcpc,
repel = TRUE, # Avoid label overlapping
show.clust.cent = TRUE, # Show cluster centers
palette = "jco", # Color palette see ?ggpubr::ggpar
ggtheme = theme_minimal(),
main = "Factor map"
)## $`1`
## v.test Mean in category Overall mean sd in category Overall sd
## Ce -2.400916 14.236667 27.547727 9.9841602 27.8314601
## Cd -2.693181 1.470833 1.664091 0.2436342 0.3602229
## Pb -3.566253 5.405833 5.980000 0.5037933 0.8082135
## p.value
## Ce 0.016354097
## Cd 0.007077392
## Pb 0.000362122
##
## $`2`
## v.test Mean in category Overall mean sd in category Overall sd
## Pb 3.700797 6.764444 5.980000 0.4474399 0.8082135
## Cd 2.685257 1.917778 1.664091 0.3517821 0.3602229
## Zn -2.243528 45.891111 59.403182 16.3646560 22.9641130
## p.value
## Pb 0.000214923
## Cd 0.007247400
## Zn 0.024862798
##
## $`3`
## v.test Mean in category Overall mean sd in category Overall sd
## Ce 3.562238 126.69 27.54773 0 27.8314601
## pH 2.377132 8.00 6.53000 0 0.6183923
## p.value
## Ce 0.0003677074
## pH 0.0174478606
Honoria_cluster <- res.hcpc$data.clust
View(Honoria_cluster)
library(tidyverse)
Honoria_cluster %>%
filter(clust=="1")## [1] "sf" "data.frame"
## Este Norte Cd Pb
## Min. :518642 Min. :9031868 Min. :1.100 Min. :4.740
## 1st Qu.:519033 1st Qu.:9031913 1st Qu.:1.350 1st Qu.:5.247
## Median :519183 Median :9032036 Median :1.575 Median :5.850
## Mean :519108 Mean :9032052 Mean :1.634 Mean :5.929
## 3rd Qu.:519221 3rd Qu.:9032134 3rd Qu.:1.863 3rd Qu.:6.525
## Max. :519556 Max. :9032343 Max. :2.550 Max. :7.240
## Zn pH Ce lng
## Min. :11.80 Min. :5.410 Min. : 4.50 Min. :-74.83
## 1st Qu.:44.79 1st Qu.:6.160 1st Qu.: 9.00 1st Qu.:-74.83
## Median :58.70 Median :6.375 Median : 13.50 Median :-74.83
## Mean :58.66 Mean :6.494 Mean : 26.11 Mean :-74.83
## 3rd Qu.:71.89 3rd Qu.:6.588 3rd Qu.: 40.02 3rd Qu.:-74.83
## Max. :99.40 Max. :8.000 Max. :126.69 Max. :-74.82
## lat geometry
## Min. :-8.758 POINT :24
## 1st Qu.:-8.758 epsg:4326 : 0
## Median :-8.757 +proj=long...: 0
## Mean :-8.757
## 3rd Qu.:-8.756
## Max. :-8.754
Honoria_elementos4 <- Honoria_elementos3 %>% st_set_geometry(NULL)
coordinates(Honoria_elementos4) <- ~Este+Norte
class(Honoria_elementos4)## [1] "SpatialPointsDataFrame"
## attr(,"package")
## [1] "sp"
## Object of class SpatialPointsDataFrame
## Coordinates:
## min max
## Este 518642 519556
## Norte 9031868 9032343
## Is projected: NA
## proj4string : [NA]
## Number of points: 24
## Data attributes:
## Cd Pb Zn pH
## Min. :1.100 Min. :4.740 Min. :11.80 Min. :5.410
## 1st Qu.:1.350 1st Qu.:5.247 1st Qu.:44.79 1st Qu.:6.160
## Median :1.575 Median :5.850 Median :58.70 Median :6.375
## Mean :1.634 Mean :5.929 Mean :58.66 Mean :6.494
## 3rd Qu.:1.863 3rd Qu.:6.525 3rd Qu.:71.89 3rd Qu.:6.588
## Max. :2.550 Max. :7.240 Max. :99.40 Max. :8.000
## Ce lng lat
## Min. : 4.50 Min. :-74.83 Min. :-8.758
## 1st Qu.: 9.00 1st Qu.:-74.83 1st Qu.:-8.758
## Median : 13.50 Median :-74.83 Median :-8.757
## Mean : 26.11 Mean :-74.83 Mean :-8.757
## 3rd Qu.: 40.02 3rd Qu.:-74.83 3rd Qu.:-8.756
## Max. :126.69 Max. :-74.82 Max. :-8.754
Grilla generada por Meuse Data:
# Hacemos una grilla regular con `SpatialPolygonsDataFrame`
proj4string(Honoria_elementos4) <- CRS("+proj=utm +zone=18 +south +datum=WGS84") #definimos proyecciónHonoria_elementos4 <- "+proj=longlat +datum=WGS84 +no_defs"
grd <- makegrid(Honoria_elementos4, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos
# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
coords = grd,
proj4string = CRS(proj4string(Honoria_elementos4))
)
# Then, visualize your clipped grid which can be used for kriging
ggplot(as.data.frame(coordinates(grd_pts))) +
geom_point(aes(x, y))
gridded(grd_pts) <- TRUE
grd_pts <- as(grd_pts, "SpatialPixels")## OGR data source with driver: ESRI Shapefile
## Source: "D:\CursoDAG\DataAnalisisGeociencias\Modulos\ModuloII\Data_Suelos\shp\Cacao_shp\Limite_Cacao..shp", layer: "Limite_Cacao."
## with 1 features
## It has 1 fields
## Integer64 fields read as strings: Id
grd <- makegrid(Limite2, n = 1000) #generamos grilla
colnames(grd) <- c("x", "y") #asignamos nombres
plot(grd$x,grd$y) #visualizamos# Next, convert the grid to `SpatialPoints` and subset these points by the polygon.
grd_pts <- SpatialPoints(
coords = grd,
proj4string = CRS(proj4string(Limite2))
)
# subset all points in `grd_pts` that fall within `spdf`
grd_pts_in <- grd_pts[Limite2, ]
grd_pts_in2 <- grd_pts_in
proj4string(grd_pts_in2) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")
gridded(grd_pts_in2) <- TRUE
grd_pts_in2 <- as(grd_pts_in2, "SpatialPixels")Honoria_elementos5 <- Honoria_elementos4
proj4string(Honoria_elementos5) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")
cd.idw = idw(Cd~1, Honoria_elementos5, grd_pts_in)## [inverse distance weighted interpolation]
## [inverse distance weighted interpolation]
m <- vgm(.59, "Sph", 874, .04) #setear bien
x <- krige(log(Cd)~1, Honoria_elementos5, grd_pts_in2, model = m)## [using ordinary kriging]
El analisis de informacion se puede realizar casi completamente en el software R y Rstudio, lo que permite automatizar el proceso en caso de estudios similares.
El análisis realizado logró mostrado como un estudio exitoso en una revista indexada de nivel Q1 en Heavy Metal Pollution and Its Effects on Agriculture en MDPI, El link es el siguiente: Agronomy MDPI.